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Abstract 



Motivated by the idea of developing a "hydrodynamic" description of spatiotem- 
poral chaos, we have investigated the defect-defect correlation functions in the de- 
fect turbulence regime of the two-dimensional, anisotropic complex Ginzburg-Landau 
equation. We compare our results with the predictions of generic scale invariance. 
Using the topological nature of the defects, we prove that defect-defect correlations 
cannot decay as slowly as predicted by generic scale invariance. We also present results 
on the fluctuations of the amplitude field A. 
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1 Introduction 



Chaos is the name given to intrinsic random behavior arising in a deterministic system [|T|. 
The simplest forms of chaos occur in systems of three or more coupled ordinary differential 
equations and in discrete mappings. A great deal is known about such systems with a small 
number of degrees of freedom displaying "temporal chaos" , where the structure of the phase 
space can be analyzed in detail. Spatially extended systems with many interacting degrees of 
freedom exhibit "spatiotemporal chaos" [0], in which the chaotic fluctuations occur on spa- 
tial scales as well as in time. Spatiotemporal chaos has also been called "weak turbulence" , 
"defect chaos", and "defect turbulence". (This is to be distinguished from fully developed 
turbulence, which typically occurs for increasing driving but fixed size L, while spatiotem- 
poral chaos occurs for fixed large driving but increasing L.) There have been a number 
of approaches used in attempting to understand spatiotemporal chaos: partial differential 
equations, coupled ordinary differential equations (in which space is discretized), coupled 
map lattices (in which space and time are discretized), and cellular automata (in which 
space, time, and dependent variable are discretized). For large systems, we may hope to use 
a statistical description of the chaotic states, borrowing concepts from equilibrium statistical 
mechanics. We wish to find simple reduced descriptions, emphasizing the collective behavior 
of many chaotic degrees of freedom, analogous to the reduced long-wavelength description 
provided by thermodynamics and hydrodynamics for equilibrium systems ^ ^ ^ |], 0. 

Typically, spatiotemporal chaos appears in nonequilibrium pattern-forming systems slight- 
ly above their threshold of instability |T^. The study of spatiotemporal chaos has been 
advanced by the development of experimental systems which are precisely controlled and 
have a large aspect ratio [0, |T^, 0, |T3|, |T^, |TB|, |T^, |T^, [T^ . Such systems have large statisti- 
cally homogeneous regions relatively free from boundary effects. A key question to address 
is whether such regions can be described in terms of hydrodynamic-like theories, focusing 
on collective behaviors and long-wavelength descriptions. We wish to address an aspect of 
this question by considering the coherent structures known as topological defects (or spirals 
or vortices) in a system that exhibits spatiotemporal chaos: the complex Ginzburg-Landau 



equation 12^, ^ |2J]. This equation describes the slowly varying amplitude and 



phase in an extended system which undergoes a supercritical Hopf bifurcation to an oscillat- 
ing and spatially uniform or oscillatory and spatially periodic state. The equation exhibits 
many interesting patterns, but we will restrict our investigation to the Benjamin-Feir un- 



stable (or defect turbulent) regime |25], in which topological defects occur in the context 
of spatiotemporal chaos. Other systems, such as Rayleigh-Benard convection |[TT], |T3 
electrohydrodynamic convection in liquid crystals [|15|, |16|, |1^ ^ capillary ripples 



cardiac tissue ||2^, chemical reactions [^, and wide aperture lasers [|T9], 0, can exhibit 
similar defect-turbulent behavior. Here, we will examine the defect-defect correlation func- 
tions and relate them to the ideas of generic scale invariance, which is a theory for describing 
nonequilibrium systems with conservation laws. We will prove that the generic predictions 
cannot be correct for topological defect correlations. Finally, we will see that our numerical 
results do not agree with the predictions of generic scale invariance ^ ^ 0, 

This paper is organized as follows. In Section ^ we discuss briefly the complex Ginzburg- 
Landau equation. The numerical technique that we use to solve this equation is discussed 
in Section E. Some interesting results that don't depend solely upon the defects in the 
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system are presented in Section ^ In Section | we briefly describe the ideas behind generic 
scale invariance, and relate these ideas to the topological defects in the complex Ginzburg- 
Landau equation. Constraints on the behavior of the defects due to their topological nature 
are discussed in Section ^. The predictions of generic scale invariance, as well as that from 
the constraints, are compared with numerical results in Section |^. Finally, we present our 
conclusions in Section |8|. 



2 The Complex Ginzburg— Landau Equation 

Perturbative analyses of microscopic equations for various pattern-forming systems yield 
complex partial differential equations which go under the name "amplitude equations". 
Considered in their own right as model dynamical systems which do not necessarily describe 
any real physical system, these equations are referred to as Ginzburg-Landau models. It 
is important to note the role of these Ginzburg-Landau equations as model equations. 
Many properties of nonequilibrium systems are encountered in these equations. Many hard 
problems, such as the existence and interaction of defects and coherent structures, or the 
appearance of chaos, can be profitably addressed in the simple framework provided by these 
equations. However, these equations provide a quantitative description of real experiments 
valid only in a small region near the transition threshold for a pattern. Far from threshold, 
only the phase of the complex field survives as a slow degree of freedom, since it describes 
a symmetry of the system. The magnitude of the complex field only becomes slow near 
threshold; far away it is just one of the many fast degrees of freedom. 

We wish to consider spatiotemporal chaos and topological defects in an equation that has 
been used both as an "amphtude equation" in the sense discussed above as well as a model 
equation to study generic features of spatially extended nonlinear dynamical systems. We 
are using the equation as a model equation for spatiotemporal chaos; its applicability to real 
physical systems is limited to the weakly nonlinear regime. The complex Ginzburg-Landau 
equation we consider is given by 

dtA = (/ii + ifi2)A - (ci + ic2)\A\^A + {h, + ih,)dlA + {hy + ib2y)d^yA (1) 

where A is a complex field in two dimensions (see for example [3^). By rescaling A' = 
,f^Ae-'^^\ t' = lilt, x' = .f^x, y' = ,f^y, = ^^by = ^, and c = ^, we obtain the 
following equation (dropping now the primes): 

dtA = A-{l + ic) I + (1 + ib,)dlA + (1 + iby)dlA. (2) 

This equation can have topological defect solutions where A = (both Re [A] and Im[y4] are 
zero) [^, ^ ^ ^ ^ ^ These defects can occur either in static arrangements 

or in dynamic ones. The dynamics states can have defects moving in a regular fashion 
(for example, all defects drifting in one direction) or they can exhibit defect turbulence, 
where defects are continuously nucleated and annihilated in pairs and are moving about 
in a chaotic, non-regular fashion. We wish to focus on the latter case, the Benjamin-Feir 
turbulent instability regime p3, 031, HI], which occurs when 1 + baC < 0. In this region of 



parameter space, all spatially periodic solutions of the complex Ginzburg-Landau equation 
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are unstable. For comparison with the ideas of generic scale invariance [31, 32, 33, 3i, 35 
we will focus on the anisotropic case by. 

The topological defects come in two varieties. The type of defect depends on how the 
phase of A changes as we travel counterclockwise once around the defect. A defect with 
a phase jump of 2tt has a topological charge of +1, while one with a jump of —2tt has a 
topological charge of —1. This is analogous to the right-handed and left-handed single- 
armed spirals in Belousov-Zhabotinskii reactions Let p+(r) equal the density of +1 
defects and P-(r) equal the density of —1 defects. We can then define a "topological" order 
parameter, p(r) = p+(r) — P-(r), which is just the density of the defects weighted by their 
topological charge. This order parameter is conserved in a system with periodic boundary 
conditions: /yp(r)dr = 0, as defects can only be created or destroyed in +/— pairs. We 
focus on the order parameter p(r) as an effective coarse-grained field, which, for the purpose 
of testing the applicability of generic scale invariance, we conjecture can be described by 
some hydrodynamic equation of motion. Figure |1| is a snapshot of the simulated system, 
showing the two types of defects, and lines of Re [A] = and Im [A] =0. 

Although the time evolution of ensembles of defects is very complicated, considerable 
progress has been made in the study of the dynamics of isolated or weakly interacting 
defects. Defects are isolated singularities in the phase equation (for the phase of A), but 
are smooth solutions of the full equations. It is an attractive idea to imagine a description 
in terms of coupled dynamics of phase and defect degrees of freedom. Progress on this idea 
has been made in the case of Rayleigh-Benard convection |^, ^ and the Kuramoto- 
Sivaskinsky equation for a review of the situation in the complex Ginzburg-Landau 
equation see |]IU|, pages 920-922]. 

The complex Ginzburg-Landau equation has also been studied in one dimension, where 
it also exhibits spatiotemporal chaos. Much effort has been made in this case to apply 
ideas from low-dimensional chaos, such as Lyapunov exponents and dimension densities 
0, ^ ^ It is not clear that such ideas will be applicable to spatiotemporal chaos 

in general, or to the two-dimensional complex Ginzburg-Landau equation with topological 
defects in particular. In one-dimension, an important question, still unsettled, is whether 



phase turbulence (turbulence without strong amplitude fluctuations) exists [0]. In contrast, 



even an apparent phase turbulence regime has not been seen in the two-dimensional complex 
Ginzbug-Landau equation. The two dimensions (one versus two-dimensional) may very well 
be fundamentally different. Numerical simulations can provide insight into the relationship 
between these regimes. 



3 Numerical Methods 

Equation (^ is a continuum equation. We wish to put it in a form that allows approximate 
solution on a computer, which requires a discretization scheme, an eigenfunction expansion 
scheme, or some combination of both. We use a technique that combines elements of both 



approaches; it is a spectral method (for a general review of this see ||53| , p^). This is a special 
form of a more general technique known as the method of weighted residuals 

Our original code was obtained from P. CouUet, L. Gil, and J. Lega (see references 
prU , pT| , P^ , |23| , 0), and is now used as the basis for simulations by several groups. The 
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original code was simply the solution to the differential equation. We have optimized the 
code, speeding it up by a factor of 3; furthermore, we have added sections for analysis of 
defects and Fourier space quantities. 

To begin our discussion of the numerical solution algorithm, we rewrite equation (|^) as 

dtA = A + B^dlA + BydlA-C\A\^A (3) 

where = 1 + ib^, By = 1 + iby, and C = 1 + ic. Now, if we knew the exact solution, 
equation would be satisfied everywhere in space. However, we can only use a finite 
number of basis functions or grid points in a numerical solution method, so the solution 
we obtain is necessarily approximate. In order to make this approximation systematic, we 
define the residual 

i?(x, t)=dtA~A + C\A\''A - B^dlA - Byd^A (4) 

for the numerical solution. In order for our numerical solution to be an accurate approxi- 
mation to the full solution, this residual needs to be small. There are a number techniques 
that accomplish this and also fall under the general name "method of weighted residuals" 

To formulate the numerical solution, the field A{'x,t) is expanded 

Nb 

A(x,t) = ^A„(t)0„(x) (5) 

a=l 

where 0a (x) is some basis function. Then we demand that 

y'i?(x,t)^^(x)dx = 0, P = 1,...,Nb (6) 

where i/^pi'x.) is our weighting function. A variety of numerical techniques follow this general 
procedure. For example, if the basis set {4>a} is the same as {ip^} then the method is known 
as a Galerkin technique. The finite element method can also be written in this general form 



55, 56 



We use the "Fourier collocation" method |53, M. Note that we are always using periodic 



boundary conditions. For other types of boundaries (such as A = at the boundaries), we 
could use Fourier sine transforms and expand in a basis set of sine waves. Our approach is 
to use 

0a(x) = e^i^-^ (7) 

where ka = (f^^i, j^f^2) = {k^, ky) for ki = — 1 and /c2 = 0, iV^ — 1. This is the 

"Fourier" part. For the "collocation", the residual is forced to be zero at the lattice points 
by choosing a set of delta functions for the weight functions 

W = (5(x - x^) = 6{x - Xj^)6{y - Vj^) (8) 
where = jidx = jijr and i/j^ = j2Sy = j2jr, with ji = 0,...,N^-1 and ja = 0, 
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At this point we define our conventions for the Fourier transform. The forward transform 

is 

^ y ji=o j2=o 



and the inverse transform is 



, ^ , ^ .ilk-, n_; J2 

= E E /^....e^-^e^-'— . (10) 

fci=0 k2=0 

We calculate the Fourier transforms in the algorithm using IBM's Engineering and Scientific 
Subroutine Library (ESSL) Fast Fourier Transform (FFT) [p8| . 

We now want to substitute equations and (H) into equation (|^). We will show this 
for the term d'^A. Start with 

1 9^A(x,t)dx = - E A,Mt) (^) e'^''^e'''''-^6ix - j,5x)5{y - j25y)dxdy. (11) 

kik2 \ ^ / 

Performing the integral and recalling the definitions of 5x, Sy, and kx gives 

EA....(t)A:^e-^e^-^. (12) 

fclfc2 



Doing similar integrals for the other terms in -R(x, t) and then using the inverse Fourier 

transform on the expressions (which amounts to dropping the terms e '^'"^ and e as 
well as the summation) gives us the equation (we write k for kik2) 



dtAkit) = (1 - Bxkl - Bykl)Ak{t) + Nk{t) (13) 

where Nk{t) is the Fourier transform of — C|y4py4 (the nonlinear part of equation (^). We 
simplify this even further by defining 

^ (1 - Bxkl - Bykl) (14) 

to give 

dtA,{t) = LkAkit) + N,it). (15) 
Multiplying both sides of equation (|1^) by an integrating factor e"^** and simplifying gives 

[e-^'=*Afc(t)] =e-^'=*iVfc(t). (16) 

We then integrate from t to t + 5t to obtain 

rt+St 

^-L,it+st)^^^^ + 6t) - e-^'=Mfe(t) = e-'^'' Nk{t')dt' (17) 

or 

rt+St 

Ak{t + 6t) = e'^'^^'Akit) + e^'=(*+'^*) e"'''^' Nk{t')dt' . (18) 
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Now, consider the integral in equation ([T8|). We approximate it as follows: 



rt+St 



-LU' 



Nkit')dt'^Nkit + 6t) 



t+5t 



-Lkt' 



dt' = Nkit + 6t) 



Q—Lkt g— Lfe{t+5i) 



.L,St 



1 



Nk{t + 5t). 



Substituting this into equation (|T8D yields 

Ak{t + 5t)=e'^^''Au{t) + 

Finally, we use the Adams-Bashforth second-order time step for A^a:(^ + ^i)'- 

5t). 

This gives us the full time step equation 



Nk{t + 5t) = ^Nk{t) - ^Nkit 



^N,{t) - ^N,{t - St) 



(19) 



(20) 



(21) 



(22) 



The solution algorithm for equation ( p2D then consists of the following steps. 

1. Calculate — CjApA in real (jij2) space. 

2. Transform the result of step 1 to Fourier (fci/c2) space to obtain Nk{t). 
3. 



Perform the time step given by equation (p2D to obtain Ak{t + 6t). Save this result to 
use in calculating y4fc(t + 26t). 

4. Invert Ak(t + 6t) to obtain y4(x, t + 6t). 

5. Save Nk{t) to use in the Adams-Bashforth part of equation (|2^ ) for calculating Ak{t + 
26t). 

Analysis of various A;-space quantities can be performed between steps 3 and 4 in the algo- 
rithm. For the initial time step we use Nk{—6t) = Nk{0). We also have started our systems 
with two different types of initial conditions. The results do not depend upon which we 
choose. One initial condition is to seed the system with two oppositely charged topological 
defects. The other is to begin with random initial conditions (for instance a uniform distri- 
bution between —0.1 and 0.1 for the field A). After starting, we have to wait for the system 
to "equilibrate" to a steady-state condition (meaning fluctuations about an average number 
of defects). This takes, for our typical parameters, on the order of 5, 000 time steps with a 
step size of 6t = 0.02. The steady-state condition is a statistically stationary state, where 
averages depend only on differences of space and time coordinates. 

By using a parallel (distributed) FFT, larger systems can be studied In the parallel 



Np 



n - 



algorithm, we distribute the data for the system over strips {i.e. the data for ji 
1), and j2 = 1, Ny exists on processor n, where Np is the total number of processors). 
The ^ FFT's for the y (j2) direction are calculated locally on each processor. We then 
perform a global transpose of the data for A by using collective communication (all to all) 
message passing. Once we have done this, the FFT's for the x (ji) direction are calculated 
locally on each processor. Then we have the result in Fourier space. To invert the FFT, we 
just reverse the process. 
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3.1 Improving the Time Step 



The result given for equation (^) represents a somewhat uncontrolled approximation that 
was made in the original code in order to minimize storage requirements. There is a more 
accurate time step (which takes somewhat more storage), which we explain in this section 
(see for a description of this process). 

We wish to find an approximation to the equation: 



(23) 



where tp = p6t = t, tp-j = tp — j6t, and tp+i = tp + I6t. We then approximate Nk{t') using 
an interpolating polynomial 



Q 

i=0 



where 



L^it')= n 



t t/f, — 



p-l 



is the interpolating polynomial. Substituting this into ( p3|) gives 



I = e^^(*+^*) E N.itp.,) f'"' e-^^*'L,(t')dt'. 

i=0 "^^p-i 



(24) 



(25) 



(26) 



For comparison with the standard Adams-Bashforth formula, we are interested in the case 
j = 0, / = 1, and q = 1. Other choices for j, /, and q give other time-stepping schemes 
{e.g. Crank-Nicholson is given by j = 1, / = 0, and q = 1). For our case the interpolating 
polynomials are 

t' - tp_i t' -t + 6t 



and 

Using these results we obtain: 



Loit') = 
Li{t') 



tp tp—i 



t' -U 



tp—i tp 



6t 

t'-t 
-6t ' 



and 



e"^^* Lo{t')dt' 



'''^\-L.t'L,{t')dt' = -- 



-2 



LkSt 



+ 1 



LkM 



Lk 





1 - 


[ V Lk5t) 


LkSt. 



(27) 
(28) 

(29) 
(30) 



The full integral (^) becomes 



Lk 



Nk{t) 
-Nkit - 5t) 



,LkSt 



1 + 



Li,6t 



1 

Lk 



+ 



LkStJ ' LUt 



(31) 
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In the limit 0, equation ( ]29|) goes to |5t and equation (|30[) goes to —^St, which 

are just the results for the standard Adams-Bashforth second-order formula. Also, when 
Lfc is held fixed and 5t — >■ 0, the leading order terms become ^St and ^St, with corrections 
of 0{6t'^), so that the original time step (equation (p2|)) is still "correct" to 0{6t). However, 
the improved time step developed in this section allows simulations to be run with a larger 
time increment 6t for the same desired numerical accuracy. This time step method has been 
used for simulations of Rayleigh-Benard convection ^T|]. Our results were obtained with the 
unimproved time step, but we present the improved version here so that others may take 
advantage of it. 



4 Results for A(x, t) 

We present some results for the field A(x, t). First, Figure ^ shows that the field-field 
correlation function decays exponentially (or perhaps faster), and shows no hints of long- 
range order (such as a power law decay might suggest). We also note that the x and y 
directions show statistically different behavior. The y direction appears to show a domain 
type structure (as in, for example, spinodal decomposition [^) for r ^ 10, because of the 
significant amount of negative correlations in that region. The x direction does not appear 
to show this behavior, as it decays smoothly to zero from its value at r = 0. 

The exponential decay of the A~A correlation function has been observed in previous 
simulations and it was suggested by these authors that the defects are responsible 

for the strong decorrelation of the field. We can merely note here that the field A(x, t) is 
not conserved in time. This lack of a conservation law would imply, within the framework of 
generic scale invariance discussed below in Section ^, that the field-field correlation function 
should decay exponentially. 

In Figure || we show the averaged power spectrum {S{kx)) = Re{A*{kx)A{0)); results are 
similar for {S{ky)). This falloff is exponential (or perhaps faster), which is to be expected 
because of the fast decay of the field-field correlation function presented in Figure |. The 
fact that the power spectrum decays rapidly to zero means that we don't have problems 
with aliasing of the Fourier transform in the numerical method. 

The probability of a particular value of Re[y4(x)] versus Re[74(x)] is shown in Figure ^. 
This figure includes both spatial averaging at a particular time, as well as time averaging 
at a particular point in order to test the system's ergodicity. The plot shows no special 
structure, and is not expected to, because principles that might give a special form, such 
as the Central Limit Theorem, do not apply to the measured quantity. Similar results are 
obtained for the probability of Im[A(x)]. We note that the two methods of averaging agree 
fairly well. We would expect that for larger systems and longer times the two results should 
agree. In Figure | we present the probability of a particular value of |A(a;)| versus |A(a;)|, 
again with both types of averaging. Values of \A\ greater than 1.0 represent shock waves 
in the system, and the probability of larger and larger values of 1^41 should decrease above 
\A\ = 1.0 When |y4| = 0, there is a defect in the system. This is expected to be rarer 
than other nonzero values, because it requires both KeA = and IraA = 0. Both of these 
expected results are seen in the data. 

In Figure we show the result for the probability of A{k) versus A{k). Non-Gaussian 
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fluctuations in Fourier variables have been measured in capillary wave fields, and this has 
been used to question the applicability of a thermodynamic description for spatiotemporal 
chaos In contrast to the results for A{x), we expect the Central Limit Theorem to 

hold for A{k). This would then predict that the probability distribution function would be 
a Gaussian. In Figure || we compare our results to a Gaussian, and we see that we certainly 
have Gaussian fluctuations. Similar results hold for other wavevectors k for which we have 
measured this quantity. We also measured this probability distribution for these fluctuations 
for a different set of parameters (c = —0.5, = hy = 40.0) and started with a small system. 
The small system showed non-Gaussian fluctuations. When we went to larger systems, the 
fluctuations became Gaussian. This just means that for the small system there were not 
enough underlying degrees of freedom for the central limit theorem to be valid when we 
averaged in the system. Note that this small system also showed quasiperiodic behavior in 
the time series of y4(x, t) for fixed x. This system had = Ly = 60. When we increased 
the size of the system to = Ly = 120, this quasiperiodic behavior disappeared and was 
replaced by the chaotic behavior seen in p4|. This larger system was the one that also 



showed the Gaussian fluctuations. 

Time averages of spatiotemporal chaos in experimental systems result in periodic spatial 
structures We also tried averaging our amplitude A(x, t) over various times, and we 

did not see any particular persistent structures emerging. These results were not, however, 
from extensive tests. It is likely that the experimental results are due to the fixed boundary 
in these systems, which constrains the system to fluctuate about certain well-defined normal 
modes. 

All of the above results for the field itself were not encouraging from the standpoint of 
illuminating special features of spatiotemporal chaos. We wish to examine another possible 
approach for describing the system, namely, generic scale invariance. 



5 Generic Scale Invariance 

In equilibrium systems, spatial correlations typically decay exponentially. For nonequilibrium 
systems, such as those with an external driving force, the situation can be quite different. 
For a nonequilibrium system with a conservation law and external noise, spatial correlation 
functions can decay algebraically. It has been suggested that this algebraic decay is expected 
to occur for a broad range of conditions, and this has been called "generic scale invariance" 
1 3l|, ^ . Some extended deterministic chaotic systems also exhibit algebraic decay 



|57| , In at least one of these examples, the chaotic fluctuations appear to play the 



same role as stochastic noise The complex Ginzburg-Landau equation would seem to 
satisfy the criterion for generic scale invariance; it shows nonequilibrium behavior since it 
cannot be derived from an underlying potential {i.e. it is non-relaxational), and in a sys- 
tem with periodic boundary conditions, the topological order parameter p(r) is conserved. 
Finally, we conjecture that the chaotic noise in our system plays the role of stochastic noise, 
again for the purpose of testing the applicability of generic scale invariance. 

With this set of conditions, we could have a hydrodynamic equation for the conserved 
order parameter (coarse grained over distances larger than the typical spacing between 
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defects) : 

dtp{r,t) = T{pir,t)} + riir,t) (32) 

where F is a general conserving operator on p, such as FqV^ + Fi(V^)^ + ^2x9^ + ^2ydy. It 
can also contain nonlinear terms {e. g. V ■ [(V^p)(Vp)] ) . The stochastic noise term rj is 
determined by: 

{Vir,t)) = (33) 

{r^{r,t)7^{r',t')) = D6{r~r')6{t-t'), (34) 

where D must be composed of differential operators for our strictly conserved order pa- 
rameter. This conserving noise term represents the effect of the chaotic fluctuations in the 
complex Ginzburg-Landau equation. There is evidence from the mapping of the Kuramoto- 
Sivashinsky equation to the Kardar-Parisi-Zhang equation |70|, |7l| and from coupled 



map lattices that this identification of spatiotemporal chaotic fluctuations with stochas- 
tic noise is not unreasonable. Note also that we are assuming that we can write down a 
local equation of motion for the order parameter p. 

For systems with nonconserving noise {i.e. D is a constant), equation is expected to 
always give rise to power law decays in the two point correlation function Gp{r) = (p(r)p(O)), 
as well as in higher order correlation functions. For nonlinear systems with conserving noise 
{e. g. D = DiV^) the situation is somewhat more complicated [0. If the system is isotropic. 



then one obtains exponential decays in Gp{r), but power law decays occur in higher order 
correlation functions. Systems which break isotropy give rise to algebraic decay in Gp{r). For 
systems with cubic symmetry, one expects Gp{r) ~ large r. For systems which 

break cubic symmetry one expects Gp{r) ~ l/r'^. We will study the last regime, a two- 
dimensional system with broken square symmetry, where generic scale invariance predicts 
that 

™(r) = (p(r)p(O)) ~ ly (35) 



for large r |]32[. We shall compare this prediction to the results from our numerics. 

The ideas behind generic scale invariance depend upon showing that nonlinearities in the 
equations are irrelevant in a renormalization group sense. This means that as the scale on 
which the system is examined grows, the importance of the nonlinearities in the dynamics 
at the larger scale is less important. This analysis is usually carried out in a perturbative 
manner, where the strength of the nonlinear term is the small expansion parameter. It 
has proven notoriously difficult to treat topological defects in a perturbative manner. An 



example of this is the Kosterlitz-Thouless transition 

We also note that generic scale invariance requires short-ranged interactions. There is 
some evidence that this occurs for the defects in the complex Ginzburg-Landau equation W^, 



74| , but collective effects might be important. Finally, the mapping of chaotic fluctuations to 
the stochastic noise could break down. In the next section, we will prove that the predictions 
of generic scale invariance cannot apply to topological defects, and therefore at least one of 
the requirements for it to be present is lacking in our system. 
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6 Topological Constraints 



We define the excess order parameter in a region to be 



5pL = 



/ (p(r + ro) - Po) dr 



(36) 



where B{L) represents a circle of radius L about a point ro, which we take as ro = due to 
translational invariance, and where po is the average order parameter (in our case 0). For 
nontopological objects, the constraint is given by Spi < aiL'^, where ai is some numerical 
constant. The excess of a nontopological object in a particular region must scale as the area 
of that region, since each individual object occupies a fixed area. For topological objects 
of the type we are studying, this constraint is different, i. e. 6pL < 0'2L, where 02 is again 
some numerical constant. The constraint arises from the fact that any excess of topological 
defects in a region must be detectable simply by traversing the perimeter of that region. 
Each topological defect has an "arm" with characteristic width that must pass through the 
perimeter of the region. Examples of this are the spiral arms of the defects in Rayleigh- 
Benard convection, extra rows of atoms for dislocations in crystals, and in our case lines of 
Re [A] = and Im [A] = 0. When a region contains the maximum excess number of defects 
allowed, each of these lines takes up a fixed amount of the perimeter of the region. Since the 
maximum excess number of topological objects scales linearly with the number of lines, and 
the number of lines scales as the perimeter of the region, we must have that the maximum 
excess number of topological objects scales as the linear size L of the region. 

If we assume that the correlation function Gp{r) decays at the same asymptotic rate 
independent of the direction of r, i.e. Gp{r) ~ f{d)g{r) where g{r) ~ for large r, 

then with this constraint we can show for two dimensions that a must be greater than 2, 
which contradicts the prediction of generic scale invariance. This result also requires that 
Jq'^ f{9)d6 7^ 0, which we expect to be true except for special cases; we discuss this point 
further in Section ^TT] below. A correlation function satisfying both assumptions occurs, for 
example, in studies of non-equilibrium conservative anisotropic lattice gases [^, |7^. To 
show that a > 2, we use the inequahty 6pL < aL. Squaring this relation yields 

/ / drdr'p(r,t)p(r',t) < a^L^. (37) 

JreB(L) Jr'eB{L) 

Now we average over the noise (or over space and time, depending upon how you wish to 
look at it) 

(p(r,t)p(r',t))^G,(r-r'). (38) 
This then gives the constraint equation 

{6pl) < a'L' (39) 

We now make the assumption that we can write Gp{r) in the form f{0)g{r) where 6 is the 
angle for r and r = |r|. Strictly speaking, this assumption is only necessary for large r. 
Given that Gp takes this form, we can write 

/ / drdr'G,(r - r') = — / /(^)dW g{R)w{R)dR, (40) 



12 



where w{R) is given by 



w 



iR)= f f S{R)6{r - r' - R)drdr' 

JrGB(L) Jr'eB{L) 



(41) 



with S{R) = 2tiR the surface area of the circle of radius R. We further will assume that 
So"" /(6')d6' 7^ 0, and now calculate a closed form expression for w{R). To do this we use the 
definition of 5(x): 



5(x). 

Plugging this into equation (^Tj) gives 

27iR 



1 



dk e 



ikx 



w(R) = / dk e-*-^ / dr e*" 

{2n)'^ J-OD JreB(L) Jr'&B{L) 



dr e-*"'. 



(42) 



(43) 



The integrals over B{L) are given by 

dr e* "" = 

reB{L) 



rdrd^e'*^"™^^ = -^Ji{kL), 



^0 



k 



(44) 



where Ji is the first order Bessel function. Substituting this result in equation (^31), we 
obtain 



w{R) = 2nRL' dke-*^'^^^^ = 4n'RL' H ^Jo ( ) Jl{k). 

3 k^ Jo k \ L 



(45) 



Performing this last integral |77] gives 



w 



(R) = inRL^ 



cos 



(46) 



We can check this by noting that Jq^ w{R)dR — vr^L^, which is just (^Jj.^b{l) 

Next, we apply w{R) to the problem at hand. Suppose that g{R) in equation ( ^OD 
exhibits its asymptotic behavior outside of some R = rmin, g{R) ~ for R > Tmin and 
also that L ^ rmin- Splitting the integral into two parts yields 



2L 







g{R)w{R)dR 



2L 



+ I g(R)w(R)dR. 
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(47) 



We consider each of these terms separately. The second term can be bounded in the following 
manner: 







g{R)w{R)dR 



< \9n 







w{R)dR. 



Expanding the integral of w{R) in descending powers of L gives 



g{R)w{R)dR 



^22 



min ijmax 



0(L). 



(48) 



(49) 
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We calculate the first integral in equation (|i7|) next. Assuming that g{R) = QqR we 
substitute equation (^61) for w{R) and write the integral as 

47rL2(2L)2-°^o C W'"" cos^^ - rf-'^il - rff^] dr/. (50) 

We again expand this integral in descending powers of L. There are no singularities due to 
the upper limit of 1, since the integrand is analytic there for all a. However, at the lower 
limit, as L oo, there can be divergences, depending on the value of a. The integrand, 
upon expansion about = 0, becomes 



(51) 



Now, suppose that a < 2. Then the integrand has an integrable singularity at = 0, and 
as L — >■ cxo we obtain some constant plus correction terms that die away for large L. For 
a = 2 we get a logarithmic divergence, and for a > 2 a diverging term L"~^ which exactly 
cancels the L^~" in front of equation (|50|) . 
The leading order behavior of {5p\) is then 



{6pl)^L'\og{L) 

{6pl) ~ (52) 



a<2 
a = 2 
a>2 

Recall that the result for generic scale invariance is the case a = 2 (see equation (|35D ). 
However, this case violates the constraint given by equation (PU]). Within our assumptions, 
this result means that for topological objects the results of generic scale invariance cannot 
hold. In fact, what we have provided is a bound on a. For topological objects, a must 
be strictly greater than 2. Generic scale invariance predicts a = 2. The simple geometric 
nature of topological objects prevents them from having correlation functions that decay 
as certain power laws. (A decay of l/(r^(logr)'^) for /5 > 1 satisfies our bounds on (5p|), 
but contradicts the predictions of generic scale invariance.) An added conclusion from the 
consideration of the topological constraints is that if the topological objects form ordered 
states, they must be of the antiferromagnetic variety {e.g. alternating + and — vortices) in 
at least one direction, in order to satisfy the topological constraint. An example of such a 
state in the complex Ginzburg-Landau equation has been seen, with defects ordering along 
chains where the defects in one chain are of opposite sign from the defects in neighboring 
chains ||f8| . 

In our analysis we have only considered the largest possible fluctuations. We expect 
these fluctuations to be rare, and hence expect a faster decay than the bound we provide. 
As an analogy, for nontopological objects the analysis presented here would predict that the 
correlation function can be at most a constant for large r; in practice, systems like spins or 
atoms have connected correlation functions that decay to zero, either as power laws or as 
exponentials. However, it is often true that inequalities in physics are saturated, especially 
in relationships between critical exponents for phase transitions. This could be the case 
here, but numerically it appears that we do not saturate this bound, as will be shown in 
Section |^ below. 
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6.1 Splitting the Correlation Function 

One further point that we should consider is the sphtting of the correlation function into 
angular and radial components. We mentioned above that this is important for our analysis. 
We cannot explicitly show this for the topological defects, but we can demonstrate this 
splitting in the context of generic scale invariance. Starting with the general result for the 
correlation function 

G{r) = J dkc^-^-^^^^ 

we then define 



G(r) = J dke"^-^ (53) 
G{r) = ^ J^'^ G{r)de (54) 



where r = (r cos 6, r sin 6) . Then 



271 Jo J r(k) 
Performing the 6 integral yields 

G{r) = J dk^Mr^kl + k^), (56) 

where Jo is the zeroth order Bessel function. Note that this expression is not in general zero. 

We can give a concrete expression, in a particular limit, for G{r). We begin with specific 
expressions for -D(k) and r(k), and ignore questions of convergence of the integrals for large 
k. Consider 

'^(■■)-£/I*-V'---'(fi|j). (57) 

The kx integral can be performed as a contour integral; the expression has poles at k^ = 
■^i^ d/cky = ±i^ky. We assume for ease of calculation the x > and y > 0; the results 
are the same for other cases. When ky > 0, we close the contour in the upper half of the 
complex kx plane, enclosing the pole at i'jky. For ky < we enclose the pole at —i'jky. The 
integral then becomes 2ni times the residue at k^ = i'yky. Then we obtain 



G(r) = — (6 - 07^) ( r dky kye-^y^'^'-'y^ + f dky kye-^y^-^'^-'A . (58) 

C7 \J0 J-oo J 



By performing the ky integrals we obtain (using x = r cos 9 and y = r sin i 



, , 2Ti{hc — ad) f d cos^ 9 — c sin^ 9 \ 1 , . 

^""^ " 7^ [{d cos^9 + c sm^9y J ^ ^ ^ 



Note that our results are for the case c ^ d, so we have for this simple example the result 
that the angular average is not zero. Also, equation ( p^D is explicitly split into a radial and 
angular component. 
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7 Numerical Results for the Defects 



We have numerically solved equation @ with periodic boundary conditions in the turbulent 
regime using a Fourier collocation code. We choose = 240, Ly = 240, = 360, and 
A^j^ = 360 (see Section ^ above). We use the parameter values c = 1.5, = —0.75, and 
by = —3.0. The time step used was 6t = 0.02. We initially "equilibrate" from a state with 
two oppositely charged defects to a state with fluctuations about some average number of 
defects. This typically takes 5000 time steps. We note that the defects do not form bound 
pairs. When a pair is created, the defects tend to move apart, and when they eventually 
annihilate, they usually do so with a defect other than their initial partner. This illustrates 
that the system is not in a Kosterlitz-Thouless bound pair phase |7^. Figure |1| shows a 
snapshot of part of our system. 

Topological defects have been studied in systems undergoing phase ordering. In these 
systems (an example is the XY model) the existence of an underlying Hamiltonian allows 
analytic progress to be made [^, ^ in the form of perturbation expansions. Some numer- 
ical work has also been performed on these systems For our system, the numerical 
results are key as we cannot form perturbation expansions due to the lack of an underly- 
ing Hamiltonian. These numerical calculations must begin with the determination of the 
location of the defects. 

To find the defects in our system, we examine the change in the phase of A as each 
plaquette (or square unit cell) on our lattice in real space is traversed counterclockwise. To 
do this, we examine the phase (/^(x, t) = tan~^ ( r^ A(x't) ) S° around a plaquette on the 

lattice and sum up the phase differences between the four points in the plaquette. A change 
of signifies that the plaquette does not contain a defect, while changes of ±27i reveal that 
a defect exists in the plaquette. We can actually speed up the defect finding process a bit. 
If all four sites in the plaquette have the same sign of Re [A] or of Im [A] , then the plaquette 
cannot contain a defect. More complicated situations can also be handled, such as where 
one corner of the plaquette is of a different sign from the other three corners for Re [A], 
while the opposite corner from the original one is of a different sign from the other three 
corners for Im [A]. In this case, the plaquette cannot contain a defect. These considerations 
avoid the costly calculation of tan~^, and can be used for a large fraction of the plaquettes 
in the lattice. 

Once we have found the defects, we can calculate n(t), the total number of defects in 
the system, and Gp{r). To do the averaging, we have run for 750,000 time steps. We only 
sample Gp(r) and n{t) every 10 time steps, because adjacent time steps are not statistically 
independent. We have calculated that (n(t)n(O)) — (n)^ ~ e~*/^, with r ~ 115 time steps. It 
has been predicted that the probability of finding a particular value of n in the system is 
given by P{n) ~ e"'-""^"^-' /^^"^ We have calculated the various moments of our distribution 
P(n), and we find (n) = 422.8 ± 0.3, cx^ = 397 ± 30, as well as a skewness of 0.014 and a 
kurtosis of —0.026, which is in good agreement with the predictions from reference [^. In 
Figure 0, we present a plot of the results for P{n) together with the fitted exponential. The 
agreement is quite remarkable. 

In Figure ^ we present the results for Gp{r) with r in both the x and y directions. For 
both directions the typical nearest neighbor is of the opposite sign: the charges are thus 



screened. Similar behavior for vortices in random wave fields has been observed 82 . In 
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Figure ^ we show a linear-log plot of \Gp{r)\. In Figure |TO]we show log-log plots of |Gp(r)|. 
We also show lines that represent the slope |Gp(r')| should have if it decayed like 
We note that at the right edge of the figure, we have reached the point where our data is 
dominated by statistical noise. It is clear that neither direction shows the expected 
decay. Our results are at variance with the predictions of generic scale invariance, as they 
must be. As was shown in the last section, the theory is not applicable to systems which 
have strong constraints placed on them due to the topological nature of the order parameter. 

The suggestion has been made that in order to see the decay the system needs 
to break x — > —x symmetry. We systematically tested this idea by adding various terms to 
the original equation (0). 

First, we added a term dd^A to the equation. This term does break the x — > —x 
symmetry. The result for the correlation function Gp{r) is given in Figure [T^. The results 
are similar to the ones seen in Figure |l^ for the results of the original equation. Note that 
the term we have added here actually only induces a drift in the entire pattern, and can 
be scaled away, so it is not surprising that it does not show behavior that differs from the 
original equation. 

Second, we added the term dd^A to the original equation. The result for Gp{r) is shown 
in Figure O. This term also breaks the x —>■ —x symmetry. However, a term of this form 



is expected to show, according to the ideas of generic scale invariance ||83|, a decay that is 
actually faster than so our results are not very informative. We cannot distinguish 

between power laws of and our results, so we must proceed to add a different term to 
the original equation. 

The third new term we added both breaks the x —x symmetry and is expected to 
show a decay, according to the predictions of generic scale invariance. This term is 
(ic^ajdApyl). The numerical results for Gp{r) are shown in Figure |T3|. We see from this that 
the correlation function does not decay like This is not surprising in the light of our 

results from Section ^, which apply to any two-dimensional topological defect in a complex 
field. 

Finally, we experimented with breaking the A — > —A symmetry in the original equation 
by adding the term ddxA"^ to the original equation. The results for Gp{r) are given in Figure 
|14|. These results no longer shows the strong anisotropy seen in the correlation function of 
the original unmodified equation. With the broken symmetry, the decay of the correlation 
function in the y direction occurs smoothly to zero from the negative values near r = 
without crossing zero. 

We have also constructed a coarse-grained order parameter field by defining 

Pcg{r) = J ds/i(s)p(r-s) (60) 

where h{s) is a coarse-graining weight function with / ds h{s) = 1. We will use the Gaussian 

we can also calculate Gp^^{r). The result is shown in Figure 
|15|. We can examine the leading behavior of this correlation function by noting that 

^..^(r) = (p.,(r)p,,(0)) = JJ dsds'h{s)h{s'){p{r - s)p(-s')). (61) 
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This can be rewritten as 



G,Jr) =jj dsds'^^--^e-^(^^+^")G,(r - (s - s')). (62) 
We now change variables ^ = s — s' and = s + s' to obtain 

^^..(r) = / / didrie-^'"'^\-^'"'^'G,{v - i). (63) 



We can perform the rj integral to obtain the final result: 

Gp.,(r) = ^ / die-^'"^'G,{v - i). (64) 
Using this final equation, we can explain the Gaussian decay seen in Figure O. Typically, 



the microscopic field Gp(r) will be given by Gq6{y) + 5'(r). Plugging this into equation (|64D 
gives 

= ^e-'V-' + _L j d^e-^V.,',,, _ J), (65) 



The behavior of the correlation function in Figure |T5] for small r is approximately Gaussian, 
due to the leading term in equation (^) as well as corrections to the overall decay due to 
the second term in equation (^). For larger r we see the domain structure of the coarse 
grained field pcg. We do not see any evidence of power-law decays in this coarse-grained 
correlation function. 

As a final comment on the numerical results presented in this section, note that the 
typical spacing between defects in our parameter regime is fairly large. As a result, we are 
not able to obtain adequate statistics for the correlation functions for large r, and we might 
not have reached the asymptotic regime where generic scale invariance should apply. Also, 
the results for the noisier [y) direction do not convincingly rule out behavior. However, 
the proof we have provided in Section ^ shows analytically that the topological nature of 
the defects precludes decay in the correlation function. 



8 Conclusions 

We have calculated a number of quantities of interest in the complex Ginzburg-Landau 
equation. These have included both properties that depend on the amplitude field A and 
properties that depend upon the defect order parameter p. The results for the field A 
provide information about the reliability of our numerics, as well as suggesting that the 
Central Limit Theorem holds for fc-space quantities. 

The results for the order parameter field p did not agree with the predictions of generic 



scale invariance [|3T| , |32| , |33| . This places a limit on the applicability of generic scale in- 
variance. Conversely, the results also also place a bound on the types of coarse-grained, 
statistical theories that can be used to describe spatiotemporally chaotic systems with topo- 
logical defects. 
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Figure 1: Snapshot of a 70 x 70 region. This is ~ 1/12 of the total area of the simulation. 
The solid lines are where Re [A] = 0, and the dashed hnes are where Im [A] — 0. Filled 
circles (•) are vortices with topological charge +1, and the open circles (o) have charge — 1. 
The typical distance between defects is 10. This picture is for c = 1.5, — —0.75, and 
by = -3.0. 

Figure 2: Re[(A*(r)A(0))] averaged over 12 different systems with — Ly — 240. The 
solid line is for r = rx while the dashed line is for r = ry. 

Figure 3: {S{kx)) = 'Re{A*{kx)A{0)) averaged over 500,000 time steps, with parameter 
values c = 1.5, bx = —1.5, and by = —1.5. The system has size = Ly = 60 with 90 
Fourier modes in each direction. Other Fourier modes show similar behavior. 

Figure 4: Prob(Re averaged over space at a particular time (solid line) and over time 

at a particular spatial location (dashed line). The spatial average is over 129,600 points, 
while the time average is over 200,000 time steps. This is for c = 1.5, b^ = —0.75, and 
by = -3.0. 

Figure 5: Prob(|/l(a;)|) averaged over space at a particular time (solid line) and over time 
at a particular spatial location (dashed line). The spatial average is over 129,600 points, 
while the time average is over 200,000 time steps. Note that values of \A\ > 1.0 represent 
shocks in the system, while the nonzero values at the limit 1^41 = represent the defects in 
the system. This is also for c 1.5, bx = —0.75, and by = —3.0. 

Figure 6: Prob(Re [^(/ca;)]) averaged over 500, 000 time steps, with parameter values c = 1.5, 
bx = —1.5, and by = —1.5. Also shown in the dashed line is a Gaussian with mean and 
standard deviation 0.08 for comparison. We have chosen a value oi k = The system 

has size Lx = Ly = 60 with 90 Fourier modes in each direction. 

Figure 7: Prob(n). This was obtained from a time series of length 70,000 time steps. 
This is for a system of size Lx — Ly — 240 with 360 Fourier modes in each direction. The 
parameter values are c = 1.5, b^ = —0.75, and by — —3.0. The dashed line is the exponential 
with mean /j, = 422.8 and = 397. 

Figure 8: Gp{r) versus r. The solid line is for the x direction, while the dashed line is for 
the y direction. Also shown is a line for G = 0. Note that G attains its asymptotic limit of 
from different sides of this hue. The parameters are c = 1.5, bx — —0.75, and by — —3.0. 

Figure 9: |Gp(r)| versus r (linear-log plot). The parameters are c = 1.5, bx — —0.75, and 
by = -3.0. 
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Figure 10: Log-log plot of |Gp(r)| versus r. The solid line corresponds to the x direction 
and the dashed line to the y direction. Also shown is a line with slope that would correspond 
to \Gp{r)\ ~ The parameters are c = 1.5, bx — —0.75, and by — —3.0. 



Figure 11: Log-log plot of |Gp(r)| versus r for the case where we have added a term 
ddxA{x,t) to the complex Ginzburg-Landau equation. Here we have d — 0.1. The solid 
line corresponds to the x direction and the dashed line to the y direction. Also shown is 
a line with slope that would correspond to |Gp(r)| ~ 1/r^. The parameters are c = 1.5, 
bx — —0.75, and by ~ —3.0. 



Figure 12: Log-log plot of |Gp(r)| versus r for the case where we have added a term 
dd^A(x,t) to the complex Ginzburg-Landau equation. Here we have d — 0.1. The solid 
line corresponds to the x direction and the dashed line to the y direction. Also shown is 
a line with slope that would correspond to |Gp(r)| ~ 1/r^. The parameters are c = 1.5, 
bx = —0.75, and by = —3.0. 



Figure 13: Log-log plot of |Gp(r)| versus r for the case where we have added a term 
ddx\A{x,t)\'^A{x,t) to the complex Ginzburg-Landau equation. Here we have d = 1.0. The 
solid line corresponds to the x direction and the dashed line to the y direction. Also shown 
is a line with slope that would correspond to |Gp(r)| ~ 1/r^. The parameters are c = 1.5, 
bx = -0.75, and by = -3.0. 



Figure 14: Log-log plot of |G'p(r)| versus r for the case where we have added a term 
ddlA{x,t)'^ to the complex Ginzburg-Landau equation. Here wc have d = 1.0. The solid 
line corresponds to the x direction and the dashed line to the y direction. Also shown is 
a hue with slope that would correspond to \Gp{r)\ ~ 1/r^. The parameters are c = 1.5, 
bx = -0.75, and by = -3.0. 



Figure 15: Gorrelation function of the coarse-grained topological order parameter. The 
coarse graining scale is cr = 10. The parameters are c— 1.5, bx = —0.75, and by — —3.0. 
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